# "Gender Identity as a Political Cue: Voter Responses to Transgender Candidates"
# Philip Edward Jones and Paul R. Brewer
# 07.19.18
# Using R version 3.3.1

library(foreign)
library(car)
library(gtools)
library(survey)
library(texreg)
library(ggplot2)
library(cowplot)
library(mediation)

library(Zelig)
library(ZeligChoice)
library(zeligverse)



#########################################################################################
# LOAD DATA AND WEIGHT SURVEYS
#########################################################################################

Study1 <- read.csv("JonesBrewer_Study1.csv")
Study2 <- read.csv("JonesBrewer_Study2.csv")

Study1W <- svydesign(ids=~1, data=Study1, weight=Study1$weight)
Study2W <- svydesign(ids=~1, data=Study2, weight=Study2$weight)




#########################################################################################
# ATE ON CANDIDATE TRAITS AND VOTE INTENTION (STUDY 1)
#########################################################################################

# IDEOLOGY
svyby(~pool.ideo, ~expt, Study1W, svymean, na.rm=T, vartype="ci")
svyttest(pool.ideo~expt, Study1W)

# REPS YOU
svyby(~pool.repsyou, ~expt, Study1W, svymean, na.rm=T, vartype="ci")
svyttest(pool.repsyou~expt, Study1W)

# OTHER TRAITS
svyttest(pool.tworthy~expt, Study1W)
svyttest(pool.moral~expt, Study1W)
svyttest(pool.authent~expt, Study1W)
svyttest(pool.focus~expt, Study1W)

# VOTE
svyby(~pool.vote, ~expt, Study1W, svymean, na.rm=T, vartype="ci")
svyttest(pool.vote~expt, Study1W)



#########################################################################################
# MEDIATION ANALYSIS FOR VOTE CHOICE (STUDY 1)
#########################################################################################

# ACME and ADE for candidate represents me
na1 <- na.omit(subset(Study1,select=c(pool.repsyou,pool.vote,expt,weight)))
med.fit <- lm(pool.repsyou~expt, data= na1, weights=na1$weight)
out.fit <- lm(pool.vote~expt+pool.repsyou, data= na1, weights=na1$weight)
med1 <- mediate(med.fit, out.fit, treat="expt", mediator="pool.repsyou", sims=1000, boot=T)
summary(med1)

# ACME and ADE for candidate ideology
na1 <- na.omit(subset(Study1,select=c(pool.ideo,pool.vote,expt,weight)))
med.fit <- lm(pool.ideo~expt, data= na1, weights=na1$weight)
out.fit <- lm(pool.vote~expt+pool.ideo, data= na1, weights=na1$weight)
med2 <- mediate(med.fit, out.fit, treat="expt", mediator="pool.ideo", sims=1000, boot=T)
summary(med2)

rm(na1, med.fit, out.fit)


#########################################################################################
# TABLE 1: HETEROGENEOUS TREATMENT EFFECTS ON VOTE CHOICE (STUDY 1 AND 2)
#########################################################################################

# Study 1
t1_ideo <- svyglm(pool.vote~expt*ideo, Study1W)
t1_pid7 <- svyglm(pool.vote~expt*pid7, Study1W)
t1_reli <- svyglm(pool.vote~expt*religiosity, Study1W)
t1_auth <- svyglm(pool.vote~expt*auth, Study1W)

# Study 2 [interaction for pid7 within each party condition]
t1_non <- svyglm(pool.vote~expt_non*pid7, Study2W)
t1_dem <- svyglm(pool.vote~expt_dem*pid7, Study2W)
t1_rep <- svyglm(pool.vote~expt_rep*pid7, Study2W)


texreg(list(t1_ideo, t1_pid7, t1_reli, t1_auth, t1_non, t1_dem, t1_rep), single.row=T, booktabs=T, caption.above=T, dcolumn=T, stars=c(.001, .01, .05, .1), symbol="\\dagger")

round(1-t1_ideo$deviance/t1_ideo$null.deviance,2)
round(1-t1_pid7$deviance/t1_pid7$null.deviance,2)
round(1-t1_reli$deviance/t1_reli$null.deviance,2)
round(1-t1_auth$deviance/t1_auth$null.deviance,2)

round(1-t1_non$deviance/t1_non$null.deviance,2)
round(1-t1_dem$deviance/t1_dem$null.deviance,2)
round(1-t1_rep$deviance/t1_rep$null.deviance,2)



#########################################################################################
# FIGURE 1: HETEROGENEOUS TREATMENT EFFECTS (STUDY 1 AND 2)
#########################################################################################

m1 <- to_zelig(t1_ideo)

i.scale <- seq(from=1, to=7, by=1)
t1_ideo_fd <- matrix(NA, ncol=3, nrow=length(i.scale))
for(i in 1:length(i.scale)){
	x1 <- setx(m1, ideo=i.scale[i], expt=0)
	x2 <- setx(m1, ideo=i.scale[i], expt=1)
	s1 <- sim(m1, x=x1, x1=x2)
	ans <- unlist(s1[["sim.out"]][["x1"]][["fd"]])
	t1_ideo_fd[i,] <- quantile(ans, probs=c(.5,.025,.975))
}
t1_ideo_fd <- as.data.frame(t1_ideo_fd)
names(t1_ideo_fd) <- c("est", "lo", "hi")
t1_ideo_fd$x <- rep(i.scale)
f1_ideo <- 	ggplot(t1_ideo_fd, aes(x=x, y=est, ymin=lo, ymax=hi)) + 
		geom_ribbon(alpha=.2) +
		geom_line(size=1) + 
		xlab("") + ylab("Treatment effect") + 
		scale_x_continuous(breaks=c(1,7), labels=c("Very\nliberal", "Very\nconservative"), minor_breaks=seq(from=2,to=6,by=1)) +
		scale_y_continuous(breaks=c(-1.5,-1,-.5,0,.5,1,1.5), minor_breaks=NULL, limits=c(-1.5,1.5)) +
		ggtitle("(a) Study 1: Ideology")	 +
		theme(plot.title=element_text(size=rel(1.2), face="bold", hjust=.5), 
				legend.position="none", 
				axis.text.x=element_text(size=10, colour="black"),
				axis.text.y=element_text(size=9, colour="black"),
				axis.title=element_text(size=11),
				strip.text=element_text(face="bold", size=11),
				strip.background=element_rect(fill="white"),
				panel.spacing=unit(2, "lines"))
f1_ideo




m1 <- to_zelig(t1_pid7)
p.scale <- seq(from=1, to=7, by=1)
t1_pid7_fd <- matrix(NA, ncol=3, nrow=length(p.scale))
for(i in 1:length(p.scale)){
	x1 <- setx(m1, pid7=p.scale[i], expt=0)
	x2 <- setx(m1, pid7=p.scale[i], expt=1)
	s1 <- sim(m1, x=x1, x1=x2)
	ans <- unlist(s1[["sim.out"]][["x1"]][["fd"]])
	t1_pid7_fd[i,] <- quantile(ans, probs=c(.5,.025,.975))
}
t1_pid7_fd <- as.data.frame(t1_pid7_fd)
names(t1_pid7_fd) <- c("est", "lo", "hi")
t1_pid7_fd$x <- rep(p.scale)
f1_pid7 <- 	ggplot(t1_pid7_fd, aes(x=x, y=est, ymin=lo, ymax=hi)) + 
		geom_ribbon(alpha=.2) +
		geom_line(size=1) + 
		xlab("") + ylab("") + 
		scale_x_continuous(breaks=c(1,7), labels=c("Strong\nDemocrat", "Strong\nRepublican"), minor_breaks=seq(from=2,to=6,by=1)) +
		scale_y_continuous(breaks=c(-1.5,-1,-.5,0,.5,1,1.5), minor_breaks=NULL, limits=c(-1.5,1.5)) +
		ggtitle("(b) Study 1: Party ID")	 +
		theme(plot.title=element_text(size=rel(1.2), face="bold", hjust=.5), 
				legend.position="none", 
				axis.text.x=element_text(size=10, colour="black"),
				axis.text.y=element_text(size=9, colour="black"),
				axis.title=element_text(size=11),
				strip.text=element_text(face="bold", size=11),
				strip.background=element_rect(fill="white"),
				panel.spacing=unit(2, "lines"))
f1_pid7


m1 <- to_zelig(t1_reli)
r.scale <- seq(from=1, to=4, by=1)
t1_reli_fd <- matrix(NA, ncol=3, nrow=length(r.scale))
for(i in 1:length(r.scale)){
	x1 <- setx(m1, religiosity =r.scale[i], expt=0)
	x2 <- setx(m1, religiosity =r.scale[i], expt=1)
	s1 <- sim(m1, x=x1, x1=x2)
	ans <- unlist(s1[["sim.out"]][["x1"]][["fd"]])
	t1_reli_fd[i,] <- quantile(ans, probs=c(.5,.025,.975))
}
t1_reli_fd <- as.data.frame(t1_reli_fd)
names(t1_reli_fd) <- c("est", "lo", "hi")
t1_reli_fd$x <- rep(r.scale)
f1_reli <- 	ggplot(t1_reli_fd, aes(x=x, y=est, ymin=lo, ymax=hi)) + 
		geom_ribbon(alpha=.2) +
		geom_line(size=1) + 
		xlab("") + ylab("") + 
		scale_x_continuous(breaks=c(1,4), labels=c("Not\nimportant", "Very\nimportant"), minor_breaks=c(2,3)) +
		scale_y_continuous(breaks=c(-1.5,-1,-.5,0,.5,1,1.5), minor_breaks=NULL, limits=c(-1.5,1.5)) +
		ggtitle("(c) Study 1: Religiosity")	 +
		theme(plot.title=element_text(size=rel(1.2), face="bold", hjust=.5), 
				legend.position="none", 
				axis.text.x=element_text(size=10, colour="black"),
				axis.text.y=element_text(size=9, colour="black"),
				axis.title=element_text(size=11),
				strip.text=element_text(face="bold", size=11),
				strip.background=element_rect(fill="white"),
				panel.spacing=unit(2, "lines"))
f1_reli


m1 <- to_zelig(t1_auth)
a.scale <- seq(from=0, to=4, by=.5)
t1_auth_fd <- matrix(NA, ncol=3, nrow=length(a.scale))
for(i in 1:length(a.scale)){
	x1 <- setx(m1, auth =a.scale[i], expt=0)
	x2 <- setx(m1, auth =a.scale[i], expt=1)
	s1 <- sim(m1, x=x1, x1=x2)
	ans <- unlist(s1[["sim.out"]][["x1"]][["fd"]])
	t1_auth_fd[i,] <- quantile(ans, probs=c(.5,.025,.975))
}
t1_auth_fd <- as.data.frame(t1_auth_fd)
names(t1_auth_fd) <- c("est", "lo", "hi")
t1_auth_fd$x <- rep(a.scale)
f1_auth <- 	ggplot(t1_auth_fd, aes(x=x, y=est, ymin=lo, ymax=hi)) + 
		geom_ribbon(alpha=.2) +
		geom_line(size=1) + 
		xlab("") + ylab("") + 
		scale_y_continuous(breaks=c(-1.5,-1,-.5,0,.5,1,1.5), minor_breaks=NULL, limits=c(-1.5,1.5)) +
		scale_x_continuous(breaks=c(0,4), labels=c("Lowest\n", "Highest\n"), minor_breaks=c(1,2,3)) +
		ggtitle("(d) Study 1: Authoritarianism")	 +
		theme(plot.title=element_text(size=rel(1.2), face="bold", hjust=.5), 
				legend.position="none", 
				axis.text.x=element_text(size=10, colour="black"),
				axis.text.y=element_text(size=9, colour="black"),
				axis.title=element_text(size=11),
				strip.text=element_text(face="bold", size=11),
				strip.background=element_rect(fill="white"),
				panel.spacing=unit(2, "lines"))
f1_auth





m1 <- to_zelig(t1_non)
p.scale <- seq(from=1, to=7, by=1)
t1_non_fd <- matrix(NA, ncol=3, nrow=length(p.scale))
for(i in 1:length(p.scale)){
	x1 <- setx(m1, pid7=p.scale[i], expt_non=0)
	x2 <- setx(m1, pid7=p.scale[i], expt_non=1)
	s1 <- sim(m1, x=x1, x1=x2)
	ans <- unlist(s1[["sim.out"]][["x1"]][["fd"]])
	t1_non_fd[i,] <- quantile(ans, probs=c(.5,.025,.975))
}
t1_non_fd <- as.data.frame(t1_non_fd)
names(t1_non_fd) <- c("est", "lo", "hi")
t1_non_fd$x <- rep(p.scale)
f1_non <- ggplot(t1_non_fd, aes(x=x, y=est, ymin=lo, ymax=hi)) + 
		geom_ribbon(alpha=.2) +
		geom_line(size=1) + 
		xlab("") + ylab("Treatment effect") + 
		scale_x_continuous(breaks=c(1,7), labels=c("Strong\nDemocrat", "Strong\nRepublican"), minor_breaks=seq(from=2,to=6,by=1)) +
		scale_y_continuous(breaks=c(-1.5,-1,-.5,0,.5,1,1.5), minor_breaks=NULL, limits=c(-1.5,1.5)) +
		ggtitle("(e) Study 2: Party ID\nNo party label")	 +
		theme(plot.title=element_text(size=rel(1.2), face="bold", hjust=.5), 
				legend.position="none", 
				axis.text.x=element_text(size=10, colour="black"),
				axis.text.y=element_text(size=9, colour="black"),
				axis.title=element_text(size=11),
				strip.text=element_text(face="bold", size=11),
				strip.background=element_rect(fill="white"),
				panel.spacing=unit(2, "lines"))
f1_non



m1 <- to_zelig(t1_dem)
t1_dem_fd <- matrix(NA, ncol=3, nrow=length(p.scale))
for(i in 1:length(p.scale)){
	x1 <- setx(m1, pid7=p.scale[i], expt_dem=0)
	x2 <- setx(m1, pid7=p.scale[i], expt_dem=1)
	s1 <- sim(m1, x=x1, x1=x2)
	ans <- unlist(s1[["sim.out"]][["x1"]][["fd"]])
	t1_dem_fd[i,] <- quantile(ans, probs=c(.5,.025,.975))
}
t1_dem_fd <- as.data.frame(t1_dem_fd)
names(t1_dem_fd) <- c("est", "lo", "hi")
t1_dem_fd$x <- rep(p.scale)
f1_dem <- ggplot(t1_dem_fd, aes(x=x, y=est, ymin=lo, ymax=hi)) + 
		geom_ribbon(alpha=.2) +
		geom_line(size=1) + 
		xlab("") + ylab("") + 
		scale_x_continuous(breaks=c(1,7), labels=c("Strong\nDemocrat", "Strong\nRepublican"), minor_breaks=seq(from=2,to=6,by=1)) +
		scale_y_continuous(breaks=c(-1.5,-1,-.5,0,.5,1,1.5), minor_breaks=NULL, limits=c(-1.5,1.5)) +
		ggtitle("(f) Study 2: Party ID\nDemocratic candidate")	 +
		theme(plot.title=element_text(size=rel(1.2), face="bold", hjust=.5), 
				legend.position="none", 
				axis.text.x=element_text(size=10, colour="black"),
				axis.text.y=element_text(size=9, colour="black"),
				axis.title=element_text(size=11),
				strip.text=element_text(face="bold", size=11),
				strip.background=element_rect(fill="white"),
				panel.spacing=unit(2, "lines"))
f1_dem



m1 <- to_zelig(t1_rep)
t1_rep_fd <- matrix(NA, ncol=3, nrow=length(p.scale))
for(i in 1:length(p.scale)){
	x1 <- setx(m1, pid7=p.scale[i], expt_rep=0)
	x2 <- setx(m1, pid7=p.scale[i], expt_rep=1)
	s1 <- sim(m1, x=x1, x1=x2)
	ans <- unlist(s1[["sim.out"]][["x1"]][["fd"]])
	t1_rep_fd[i,] <- quantile(ans, probs=c(.5,.025,.975))
}
t1_rep_fd <- as.data.frame(t1_rep_fd)
names(t1_rep_fd) <- c("est", "lo", "hi")
t1_rep_fd$x <- rep(p.scale)
f1_rep <- ggplot(t1_rep_fd, aes(x=x, y=est, ymin=lo, ymax=hi)) + 
		geom_ribbon(alpha=.2) +
		geom_line(size=1) + 
		xlab("") + ylab("") + 
		scale_x_continuous(breaks=c(1,7), labels=c("Strong\nDemocrat", "Strong\nRepublican"), minor_breaks=seq(from=2,to=6,by=1)) +
		scale_y_continuous(breaks=c(-1.5,-1,-.5,0,.5,1,1.5), minor_breaks=NULL, limits=c(-1.5,1.5)) +
		ggtitle("(g) Study 2: Party ID\nRepublican candidate")	 +
		theme(plot.title=element_text(size=rel(1.2), face="bold", hjust=.5), 
				legend.position="none", 
				axis.text.x=element_text(size=10, colour="black"),
				axis.text.y=element_text(size=9, colour="black"),				
				axis.title=element_text(size=11),
				strip.text=element_text(face="bold", size=11),
				strip.background=element_rect(fill="white"),
				panel.spacing=unit(2, "lines"))
f1_rep




theme_set(theme_gray())
plot_grid(f1_ideo, f1_pid7, f1_reli, f1_auth,     ncol=4, rel_heights=c(.15,1,1,1,1), scale=.93)
ggsave("f1a.tiff", width=16, height=7, units='in', dpi=800)

plot_grid(f1_non, f1_dem, f1_rep,     ncol=3,  scale=.93)
ggsave("f1b.tiff", width=12, height=7, units='in', dpi=800)

dev.off()




#########################################################################################
# ATE ON CANDIDATE TRAITS AND VOTE INTENTION (STUDY 2)
#########################################################################################

# STUDY 2 IDEOLOGY
svyby(~pool.ideo, ~expt, Study2W, svymean, na.rm=T, vartype="ci")
svyttest(pool.ideo~expt, Study2W)

# STUDY 2 REPS YOU
svyby(~pool.repsyou, ~expt, Study2W, svymean, na.rm=T, vartype="ci")
svyttest(pool.repsyou~expt, Study2W)

# STUDY 2 OTHER TRAITS
svyttest(pool.tworthy~expt, Study2W)
svyttest(pool.moral~expt, Study2W)
svyttest(pool.authent~expt, Study2W)
svyttest(pool.focus~expt, Study2W)

# STUDY 2 VOTE
svyby(~pool.vote, ~expt, Study2W, svymean, na.rm=T, vartype="ci")
svyttest(pool.vote~expt, Study2W)




# STUDY 2 IDEOLOGY by PARTY
svyby(~pool.ideo, ~expt_dem, Study2W, svymean, na.rm=T, vartype="ci")
svyttest(pool.ideo~expt_dem, Study2W)
svyby(~pool.ideo, ~expt_rep, Study2W, svymean, na.rm=T, vartype="ci")
svyttest(pool.ideo~expt_rep, Study2W)


# COMPARE TRANS R TO CIS D
Study2$expt3 <- recode(Study2$expt_conds, "'6 TreatmentRep'=1; '3 ControlDem'=0; else=NA")
Study2W <- svydesign(ids=~1, data= Study2, weight= Study2$weight)
svyby(~pool.ideo, ~expt3, Study2W, svymean, na.rm=T, vartype="ci")
svyttest(pool.ideo~expt3, Study2W)



# STUDY 2 REPS YOU by PARTY
svyby(~pool.repsyou, ~expt_dem, Study2W, svymean, na.rm=T, vartype="ci")
svyttest(pool.repsyou~expt_dem, Study2W)
svyby(~pool.repsyou, ~expt_rep, Study2W, svymean, na.rm=T, vartype="ci")
svyttest(pool.repsyou~expt_rep, Study2W)



# STUDY 2 VOTE by PARTY
svyby(~pool.vote, ~expt_dem, Study2W, svymean, na.rm=T, vartype="ci")
svyttest(pool.vote ~expt_dem, Study2W)
svyby(~ pool.vote, ~expt_rep, Study2W, svymean, na.rm=T, vartype="ci")
svyttest(pool.vote ~expt_rep, Study2W)

